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1 Introduction 

In the analysis of MEG data one often faces the prob¬ 
lem that noise from biological or technical origins 
(e.g. alpha activity or interference from the power 
line, respectively) is corrupting the measurements. 
We present a case study where we analyze the effects 
of artifact removal for a well-known experimental set¬ 
ting: measurements of somatosensory evoked fields 
(SEF, N20). We compare a classical signal process¬ 
ing approach to the recently developed independent 
component analysis (ICA) technology [9, 11]. The 
specific data set studied is an attractive testbed since 
the signal of interest (N20) is relatively strong, but 
contaminated by a 150 Hz component due to power 
line interference. 

2 Methods 

2.1 Data 

The right median nerve was stimulated electrically 
over 12000 epochs, while the magnetic field above 
the contralateral somatosensory cortex was measured 
by the Berlin 49 channel planar gradiometer system 
placed in a magnetically shielded room. A sampling 
rate of 2 kHz and an inter-stimulus interval (ISI) of 
333 msec (to avoid steady state effects) were used (see 
Fig. 1 and 2). 




Figure 2: SEF: Localization of an ECD at the time 
instance of N20 using a realistic volume conductor 
model. 


2.2 ICA model 

Due to the fact that magnetic fields of different 
bioelectric current sources superimpose linearly, the 
measured values of the SQUID-sensor array can be 
modeled as a linear combination of component vec¬ 
tors 

x(i)=As(f), (1) 

where x = [x\, ... ,x m ] T , s = [si,... , s n ] T , m > 
n. For independent component analysis we assume 
that the observed signals x(f) are linear mixtures of 
n underlying sources s(t), that are mutually statisti¬ 
cally independent, i.e. their joint probability density 
function factorizes. 

Within these assumptions one can separate the data 
x(f) into independent components u(f) = Wx(t). 
This recovers the original sources s(t) from the ob¬ 
served mixtures up to scaling and permutation. As 
both the mixing process A and the sources s(t) are 
unknown, this technique is called blind source sepa¬ 
ration [3]. 


Figure 1: SEF: averaged time courses of the magnetic 
field for all 49 channels. 



2.3 Three Algorithms 

In the following we will briefly review three represen¬ 
tative types of source separation algorithms that take 
different approaches to achieve a demixing. 

A substantial amount of research has been conducted 
on algorithms using higher-order statistics for esti¬ 
mation of ICA [3]. For off-line (batch) computa¬ 
tion, Cardoso et al. [2] developed the JADE algo¬ 
rithm based on the (joint) diagonalization of matri¬ 
ces obtained from ‘parallel slices’ of the fourth-order 
cumulant tensor. Maximizing the kurtosis of the out¬ 
put signals was also proposed by Hyvarinen and Oja 
[4]. They developed an algorithm termed FastICA 
that utilizes a fixed-point iteration to optimize a con¬ 
trast function that measures the distance of the source 
probability distributions from a Gaussian distribution 
[3, 4], 

In matrix notation FastICA takes the form 

w' = W + r[diag(-A) + £{p(u)u r }]W, 

where u = Wx, fa — E{uig(ui)} and T = 
diag(l /(fa—E{g'(ui)})) where g(u,) is a non-linear 
contrast function. If one would use g{ui) — uf the 
Kurtosis is optimized (as in JADE). Therefore we 
chose the hyperbolic tangent and the Gaussian func¬ 
tion as non-linearities g(-) in our experiments with 
FastICA. 

The two methods mentioned above utilize higher- 
order statistics to exploit the non-Gaussian distribu¬ 
tion of the sources to achieve a separation. In con¬ 
trast the TDSEP (Temporal Decorrelation SEPara- 
tion) algorithm [10, 11] relies on distinctive spec¬ 
tral/temporal characteristics of the sources using only 
second-order statistics in the form of (time-delayed) 
correlation matrices (see also [6, 1]). 

The TDSEP algorithm makes use of the property that 
the cross-correlation functions of mutually indepen¬ 
dent signals are approximately zero. Assuming fur¬ 
ther that the signals s(f) have a temporal structure 
i.e. a ’non-delta’ auto-correlation function all time- 
delayed correlation matrices i? T ( s ) should be non¬ 
zero diagonal matrices. This knowledge is used to 
calculate the unknown mixing matrix A in eq. (1) 
by a simultaneous diagonalization of a set of corre¬ 
lation matrices 1 R r ( x ) — (x(f)x T (f — t)) for differ¬ 
ent choices of t. Since the mixing model in eq. (1) is 
simply a linear transformation, we can substitute x(f) 
by As(f) and get: 

-Kt(x) = (As(f) (A s (t - t)) t ) = Ai? r(s) A r . (2) 

'here {■) denotes the time average 


For the special case of two lagged correlation matri¬ 
ces, e.g. t = 0 and r / 0 one can achieve a joint di¬ 
agonalization by solving the general eigenvalue prob¬ 
lem {R-,^o( x )R r ioi ^ A = AA [6]> 

The quality of the signal separation depends on the 
choice of t [10] because for some r the eigenvalue 
problem can be degenerated. However, solving eq. (2) 
for several 2 r by simultaneous diagonalization elimi¬ 
nates this problem. An approximate simultaneous di¬ 
agonalization of several matrices can be performed 
in two steps: (1) sphering and (2) a number of Ja¬ 
cobi rotations (orthogonal transformation). First, the 
sphering operation W — Rr= o(x) _ 2 transforms the 
covariance matrix of z(i) = Wx(f) into the iden¬ 
tity matrix. Then the remaining set of time-delayed 
correlation matrices i? T ( z ) can be diagonalized sub¬ 
sequently by a unique orthogonal transformation Q, 
since in the new basis z all degrees of freedom left 
are rotations. For details we refer to [2, 1], 
Concatenation of both transforms finally yields an es¬ 
timate of the mixing matrix A = W -1 Q, which has 
to be inverted to get the demixing matrix W = A -1 . 

2.4 Artifact reduction using ICA 

ICA has been a sucessful technique for artifact reduc¬ 
tion in EEG and MEG [5, 9, 11]. Our ICA-based 
artifact reduction procedure consist of the following 
steps: First, a sphering and compression by PCA is 
applied. For our data we reduced the 49 input dimen¬ 
sions to 15 based on the Eigenvalue spectrum shown 
in Fig. 3a. Then we decompose the transformed data 
into independent components by an ICA algorithm. 
In a next step we decide which components corre¬ 
spond to artifactual or relevant signals on the basis 
of prior knowledge. Finally we project the previously 
selected components of interest back to sensor space 
and by doing so we obtain a set of cleaned measure¬ 
ments. 

Another option is to remove the artifact fields by mak¬ 
ing use of their estimated spatial structure (contained 
in the columns of the mixing matrix A) by Signal- 
Space Projection (SSP) [7, 8]. In case of multiple 
artifacts the whole space spanned by these artifacts, 
which we will refer to as ’artifact space’, has to be 
projected out. The essential requirement for applying 
SSP is that the unwanted fields are to be known (up to 
multiplicative constants) which is fulfilled if the ICA 
merely finds the artifact space correctly. 


“In the experiments we chose r = 0..5. 




2.5 Performance Evaluation 


We make use of the fact that the ISI is relatively large 
and fit the parameters of a notch filter on the last part 
(> 200ms) of the averaged data (see Fig. 1) and ex¬ 
trapolate this to the first part. This corrected data 
is refered to as "gold standard". To assess the per¬ 
formance of the different artifact reduction methods, 
the normalized deviation from the "gold standard" is 
used. As a second criteria of interest we compare the 
goodness-of-fit of an equivalent current dipole (ECD) 
model for a realistic volume conductor (calculations 
were performed with the CURRY software package). 

3 Results 

Predominantly two artifactual 150 Hz components 
(Fig. 4) were identified by the ICA decomposition in 
a fully automatic, data-driven manner. Fig. 3b clearly 
shows that the energy present in almost all 49 chan¬ 
nels in the 150Hz band in the original data (top graph) 
is condensed to fewer ones using PCA (lower left). 
Finally in the ICA basis only two 150 Hz components 
persist (lower right). Fig. 5 reveals that ICA tech¬ 
niques eliminating these artifact fields yield si mi lar 
results - only 2-3% deviation - as the "gold standard" 
notch filtering approach, provided that the sample size 
is sufficiently large i.e. > 400 (see Fig. 5). Note that a 
crude (non-adaptive) notch filter yields a constant er¬ 
ror of approximately 8%, which is much worse than 
for the ICA based algorithms. To obtain Fig. 6, we 
project the data orthogonal to the artifact subspace 
and give the deviation from the accordingly projected 
gold standard data, which yields remarkably better 
results, especially for small sample sizes. Here we 
can obtain an improvement in deviation of a factor 10 
compared to the simple subtraction approach shown 
in Fig. 5. 

The goodness-of-fit of the dipole model is very high 
for all methods (even for the original data), due to the 
high signal to noise ratio at N20. Interestingly, if ICA 
is used, one achieves accuracies of larger than 99.8%. 

4 Discussion 

In our testbed scenario, we have access to the "ground 
truth" since we know that mainly a 150Hz distortion 
is present, and due to the large inter-stimulus inter¬ 
val we are able to adapt parameters of a notch filter 
effectively. In the general setting where no such a pri¬ 
ori information is available, it is apparent that such 



Figure 3: a) Eigenvalue spectrum of the covariance 
matrix of the data, b) relative power in the 150 Hz 
band for original data, PCA and ICA, respectively. 



Figure 4: Two components showing a strong peak at 
150 Hz were obtained by ICA (here the TDSEP al¬ 
gorithm was applied to a 15-dimensional subspace of 
the data). 



Figure 5: Deviation from the gold standard for differ¬ 
ent artifact removal methods using subtraction. 
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Figure 6: Deviation from the gold standard for 1CA- 
based methods using projection. 



Figure 7: Goodness-of-fit for an ECD using different 
methods. 


notch filtering would inevitably fail while ICA can 
still yield reasonable results. Furthermore, the com¬ 
bination of ICA with SSP to eliminate artifacts has 
proven superior to methods that merely subtract the 
noise. Our artifact reduction methods are fully auto¬ 
mated, and they could be applied to any multichannel 
data, provided that there exist criteria to spot the arti¬ 
facts. 

Finally, for the ICA-based methods, e.g. for TDSEP 
(cf. Fig. 6), shorter inter-stimulus intervals would be 
possible without losses in localization accuracy. 
Future research will furthermore focus on strategies 
for combining projection techniques and ICA incor¬ 
porating a-priori knowledge. 
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